﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.ComponentModel;
using System.Configuration;
using System.IO;
using System.Threading;
using MathNet.Numerics.Random;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.LinearAlgebra.Generic;
using MathNet.Numerics.LinearAlgebra.Single;
using System.Collections;
using System.Diagnostics;
using System.Threading.Tasks;


namespace HPSO
{
    public class HSPO
    {

        int[][] all_best = null;
        int[][] all_swarm = null;
        int[] all_particalParentIndex = null;
        Matrix[] all_velocity = null;
        long[] all_fitness = null;
        long[] all_bestfitness = null;
        Nullable<float> m_Inertia = null;
        Nullable<float> m_Cof1 = null;
        Nullable<float> m_Cof2 = null;
        Nullable<float> m_Cof3 = null;
        Nullable<float> m_MaxVelocity = null;
        Nullable<float> m_MinVelocity = null;
        Nullable<int> m_height = null;
        Nullable<int> m_degree = null;
        Nullable<long> m_bestSystemFitness = null;
        List<long> bestSystemfittness = null;
        List<int> iterationNumber = null;
        List<int[]> System_myBestResult = null;
        int[] systemBestResult = null;
        Nullable<int> m_ParticleSize = null;
        int[] m_Particlelevel = null;
        AbstractRandomNumberGenerator R = new MersenneTwister();
        BackgroundWorker backgroundWorker;
        Nullable<int> functionCounter = null;
        List<int> System_functionCounter = null;
        bool fulllocal = false;
        bool singlelocal = false;
        bool simulatedAnnealing = false;
        bool NewbiasedCondition = false;
        bool singletabou = false;
        bool Fulltabou = false;
        bool isParrallel = false;
        double temperture;
        string DirectoryName;
        
        /// <summary>
        /// this is to intialize the hierarchy best particle swarm optmization
        /// </summary>
        /// <param name="directoryName">directory name to save the results</param>
        /// <param name="filename">the problem file name</param>
        /// <param name="bw">background worker that this object will use</param>
        /// <param name="singlelocalsearch">to allow single local search</param>
        /// <param name="newbiased">to allow the biased term</param>
        /// <param name="fulllocalsearch">to allow all particles to do local search</param>
        /// <param name="parrallel">to run the particles in parallel</paraparam>
        /// <param name="tabuSearchSingle">to run tabu search on single partile</param>
        public HSPO(string directoryName,BackgroundWorker bw,string filename, bool fulllocalsearch, bool singlelocalsearch, bool newbiased, bool parrallel, bool tabuSearchSingle)
        {
            DirectoryName = directoryName;
            backgroundWorker = bw;
            m_FileNameRead = filename;
            singlelocal = singlelocalsearch;
            fulllocal = fulllocalsearch;
            NewbiasedCondition = newbiased;
            singletabou = tabuSearchSingle;
            isParrallel = parrallel;
            
        }
        /// <summary>
        /// gets the needed configuration of the system.
        /// </summary>
        private void BuildSwarmVariables()
        {
            //reads the data from the configuration file
            m_height = int.Parse(ConfigurationManager.AppSettings.Get("FirstLevelHeight"));
            m_degree = int.Parse(ConfigurationManager.AppSettings.Get("FirstLevelDegree"));
            m_Cof1 = float.Parse(ConfigurationManager.AppSettings.Get("Firstcof1"));
            m_Cof2 = float.Parse(ConfigurationManager.AppSettings.Get("Firstcof2"));
            m_Inertia = float.Parse(ConfigurationManager.AppSettings.Get("FirstInertia"));
            m_MaxVelocity = int.Parse(ConfigurationManager.AppSettings.Get("FirstMaxVelocity"));
            m_MinVelocity = int.Parse(ConfigurationManager.AppSettings.Get("FirstMinVelocity"));
            //Start to intialize the objects used in the algorithms
            m_Particlelevel = new int[m_height.Value + 1];
            m_ParticleSize = 0;
            for (int i = 0; i <= m_height; i++)
            {
                m_ParticleSize += (int)Math.Pow(m_degree.Value, i);
                m_Particlelevel[i] = (int)Math.Pow(m_degree.Value, i);
            }
            all_best = new int[m_ParticleSize.Value][];
            all_fitness = new long[m_ParticleSize.Value];
            all_bestfitness = new long[m_ParticleSize.Value];
            all_swarm = new int[m_ParticleSize.Value][];
            all_velocity = new Matrix[m_ParticleSize.Value];
            all_particalParentIndex = new int[m_ParticleSize.Value];
            functionCounter = 0;
            System_functionCounter = new List<int>();
            System_myBestResult = new List<int[]>();
        }

        string m_FileNameRead;
        string[] m_Lines;
        int m_VariablesNumber;
        int[][] m_B;
        int[][] m_A;
        Matrix m_MatrixB;
        Matrix m_MatrixA;
        Matrix m_Matrixhurstic;

        /// <summary>
        /// this function basically translate and build the problem in arrays and variables.
        /// </summary>
        public void Translate()
        {
            ReadFile();
            GetSize();
            GetArrays();
        }
        /// <summary>
        /// read the problem file.
        /// </summary>
        void ReadFile()
        {
            FileStream FS = new FileStream(m_FileNameRead, FileMode.Open, FileAccess.Read);
            StreamReader SR = new StreamReader(FS);
            string File = SR.ReadToEnd();
            //it splits the data in the file to get individual numbers
            m_Lines = File.Split(new char[] { '\n', '\r', ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
            SR.Close();
            FS.Close();
        }
        /// <summary>
        /// get the size of the problem through the first line
        /// </summary>
        void GetSize()
        {
            m_VariablesNumber = int.Parse(m_Lines[0]);
        }
        /// <summary>
        /// convert the readed data to two arrays.
        /// and create the hursitic matrix for the velocity 
        /// </summary>
        void GetArrays()
        {
            
            stopwatch_HuristicTerm.Start();
            //Intialize the objects for the problem.
            int counter = 1;
            m_B = new int[m_VariablesNumber][];
            m_A = new int[m_VariablesNumber][];
            //Intialize the objects to calculate the huristic Matrix
            m_MatrixA = new DenseMatrix(m_VariablesNumber, 1, 0);
            m_MatrixB = new DenseMatrix(1, m_VariablesNumber, 0);
            m_Matrixhurstic = new DenseMatrix(m_VariablesNumber);

            for (int i = 0; i < m_VariablesNumber; i++)
            {
                m_A[i] = new int[m_VariablesNumber];

                for (int j = 0; j < m_VariablesNumber; j++)
                {
                    m_A[i][j] = int.Parse(m_Lines[counter]);
                    m_MatrixA[i, 0] += (float)m_A[i][j];
                    counter++;
                }

            }
            for (int i = 0; i < m_VariablesNumber; i++)
            {
                m_B[i] = new int[m_VariablesNumber];
                for (int j = 0; j < m_VariablesNumber; j++)
                {
                    m_B[i][j] = int.Parse(m_Lines[counter]);
                    m_MatrixB[0, i] += (float)m_B[i][j];
                    counter++;
                }
            }

            //multiplication of both matrixes to create the hurstic variable 
            m_Matrixhurstic = (Matrix)m_MatrixA.Multiply(m_MatrixB);

            stopwatch_HuristicTerm.Stop();            
        }
        /// <summary>
        /// this function convert a result to a matrix
        /// </summary>
        /// <param name="result">the solution in hand</param>
        /// <returns>the matrix that represent it</returns>
        public Matrix ConverttoMatrix(int[] result)
        {
            Matrix Matrix = new DenseMatrix(m_VariablesNumber,m_VariablesNumber,0);
            for (int i = 0; i < m_VariablesNumber; i++)
            {
                Matrix[i, result[i]] = 1;
            }
            return Matrix;
        }
        /// <summary>
        /// generate a random solution to start the swarm with
        /// </summary>
        public int[] GenerateRandomSolution()
        {
            int[] result = new int[m_VariablesNumber];
            for (int i = 0; i < m_VariablesNumber; i++)
            {
                result[i] = i;
            }
            for (int i = 1; i < m_VariablesNumber; i++) 
            {
                result = swap(i-1,R.Next(i-1,m_VariablesNumber),result);
            }
            return result;
        }
        /// <summary>
        /// helper function to swap two indexes and returned the resulted array
        /// </summary>
        /// <param name="index">first item</param>
        /// <param name="index2">second item</param>
        /// <param name="result">result before swapping items</param>
        /// <returns>results after swapping items</returns>
        private int[] swap(int index, int index2, int[] result)
        {
            int temp = result[index];
            result[index] = result[index2];
            result[index2] = temp;
            return result;
        }
        /// <summary>
        /// calculate the fittness of the result array
        /// </summary>
        /// <param name="result">the array to calculate its fittnes</param>
        public long Total_Evaluate_For_ThreadLocal(int[] result)
        {
            long total = 0;

            for (int i = 0; i < m_VariablesNumber; i++)
            {
                for (int j = 0; j < m_VariablesNumber; j++)
                {
                    total += (m_A[i][j]) * (m_B[result[i]][result[j]]); 
                }
            }
            return total;
        }

        int iteration = 0;
        int failures;
        long intializeSystem = new long();
        long huristicTerm = new long();
        long localSearch = new long();
        long runs = new long();
        long tabuSearch = new long();
        long totalRun = new long();
        Stopwatch stopwatch_intializeSystem = new Stopwatch();
        Stopwatch stopwatch_HuristicTerm = new Stopwatch();
        Stopwatch stopwatch_LocalSearch = new Stopwatch();
        Stopwatch stopwatch_Runs = new Stopwatch();
        Stopwatch stopwatch_totalRun = new Stopwatch();
        Stopwatch stopwatch_TabuSearch = new Stopwatch();
        /// <summary>
        /// this is the run system that runs the pso system 4000 times
        /// </summary>
        public void RunSystem(object sender, DoWorkEventArgs e)
        {
            
            iteration = 0;
            //start the system intialize timer
            stopwatch_intializeSystem.Start();
            //start the total run timer
            stopwatch_totalRun.Start();
            // this function creates the object and the herarchy needed to run the PSO algorithm
            buildSystem();
            // this function create the random start point of the pso variables
            intializeSwarms();
            // Execution Time for the Instance
            TimeSpan ExecutionTime;
            // decide on a run time for the size
            if (m_VariablesNumber <= 40)
                ExecutionTime = new TimeSpan(1, 0, 0);
            else
                ExecutionTime = new TimeSpan(2, 0, 0);

            //Stop the system intiaization timer
            stopwatch_intializeSystem.Stop();
            // this is to save the variable -for the initization time
            intializeSystem = stopwatch_intializeSystem.ElapsedTicks;
            //start the run timer
            stopwatch_Runs.Start();

            failures = 10 * m_VariablesNumber;
            // run the system for 20 000 iteration
            stopwatch_Runs.Start();
            // run the system for 20 000 iteration
            while (stopwatch_totalRun.Elapsed < ExecutionTime)
            {
                failures = failures * 2;
                run();
                iteration++;
            }
            //stop both timers
            stopwatch_Runs.Stop();
            stopwatch_totalRun.Stop();
            // get the results of all timers in the system
            runs = stopwatch_Runs.ElapsedTicks;
            localSearch = stopwatch_LocalSearch.ElapsedTicks;
            huristicTerm = stopwatch_HuristicTerm.ElapsedTicks;
            totalRun = stopwatch_totalRun.ElapsedTicks;
            tabuSearch = stopwatch_TabuSearch.ElapsedTicks;
            // restart timers
            stopwatch_Runs.Restart();
            stopwatch_totalRun.Restart();
            //report produced for this run of PSO algorithm
            reportProgree();
            //clear the data
            Dispose();
        }
        /// <summary>
        /// this function is to build the arrays and varibales and itialize the arrays
        /// </summary>
        private void buildSystem()
        {
            BuildSwarmVariables();
            intializeheirarchy();
            Translate();
            intializeSwarms();
        }
        /// <summary>
        /// this is to intialize a heirarchy pointers to know each particle parent through the process
        /// </summary>
        private void intializeheirarchy()
        {
            int bestpartical = -1;
            int parentCounter = 0;
            for (int i = 0; i < m_ParticleSize; i++)
            {

                if (i == 0)
                {
                    all_particalParentIndex[i] = bestpartical;
                }
                else
                {
                    for (int l = 0; l < m_degree; l++)
                    {
                        all_particalParentIndex[i] = parentCounter;
                        i++;
                    }
                    i--;
                    parentCounter++;
                }

            }
        }
        /// <summary>
        /// this is to intialize a variables of the swarm like fittness and best solution and velocity matrix
        /// </summary>
        private void intializeSwarms()
        {
            IContinuousDistribution distribution = new Cauchy();
            m_bestSystemFitness = int.MaxValue;
            iterationNumber = new List<int>();
            bestSystemfittness = new List<long>();
            for (int i = 0; i < m_ParticleSize; i++)
            {
                int[] randomResult = GenerateRandomSolution();
                all_swarm[i] = randomResult;
                all_best[i] = randomResult;
                all_bestfitness[i] = int.MaxValue;
                all_velocity[i] = new DenseMatrix(m_VariablesNumber);
                all_velocity[i] = (Matrix)all_velocity[i].Random(m_VariablesNumber, m_VariablesNumber, distribution);
                calculateFittness(i);
            }
        }
        /// <summary>
        /// this is function to calculate fittness and set the properties like the best fittness of the particle, its best result and if it is the best
        /// fittness ever found by swarm
        /// </summary>
        /// <param name="particalIndex">particle index</param>
        private void calculateFittness(int particalIndex)
        {
            functionCounter++;
            //get the fitness of this particle
            long fitnessPartical = (long)Total_Evaluate_For_ThreadLocal(all_swarm[particalIndex]);
            //check this fittness with the best fittness of this particle
            if (fitnessPartical < all_bestfitness[particalIndex])
            {
                //if its better fittness override the best value and best fittness of this particlee
                // with the new value
                all_best[particalIndex] = all_swarm[particalIndex];
                all_bestfitness[particalIndex] = fitnessPartical;
                //check this fittness with the system best fittness of the system
                if (fitnessPartical < m_bestSystemFitness)
                {
                    //override the best fittness of the PSO algo and save the best fittness and iteration number
                    m_bestSystemFitness = fitnessPartical;
                    bestSystemfittness.Add(fitnessPartical);
                    iterationNumber.Add(iteration);
                    System_functionCounter.Add(functionCounter.Value);
                    System_myBestResult.Add(all_best[particalIndex]);
                    systemBestResult = all_best[particalIndex];
                }
            }
            //set the fittness of the particle with new value
            all_fitness[particalIndex] = fitnessPartical;
        }
        /// <summary>
        /// this is the run of the one iteration where all particles run in parallel
        /// </summary>
        public void run()
        {
            if (isParrallel)
            {
                Parallel.For
                (0, m_ParticleSize.Value, j =>
                {
                    //get veolcity
                    calculateVelocity(j);
                    // get new position through new velocity
                    reCalculate(j);
                    // if this system is full local search particle
                    // run all particles local search step
                    if (fulllocal)
                    {
                        localsearch(j);
                    }
                    // if this system is single local search particle check the index if the top particle 
                    // run local search
                    if (singlelocal)
                    {
                        if (j == 0)
                        {
                            localsearch(j);
                        }
                    }
                    // if this system is single tabu search particle check the index if the top particle 
                    // run tabu search
                    if (singletabou)
                    {
                        if (j == 0)
                        {
                            tabu_search(j, m_VariablesNumber * 8, m_VariablesNumber * m_VariablesNumber * 5, failures);
                        }
                    }
                    // recalculate the new fittness of this particle
                    calculateFittness(j);
                });
            }
            else
            {
                for (int j = 0; j < m_ParticleSize.Value; j++)
                {
                    //get veolcity
                    calculateVelocity(j);
                    // get new position through new velocity
                    reCalculate(j);
                    // if this system is full local search particle
                    // run all particles local search step
                    if (fulllocal)
                    {
                        localsearch(j);
                    }
                    // if this system is single local search particle check the index if the top particle 
                    // run local search
                    if (singlelocal)
                    {
                        if (j == 0)
                        {
                            localsearch(j);
                        }
                    }
                    // if this system is single tabu search particle check the index if the top particle 
                    // run tabu search
                    if (singletabou)
                    {
                        if (j == 0)
                        {
                            tabu_search(j, m_VariablesNumber * 8, m_VariablesNumber * m_VariablesNumber * 5, failures);
                        }
                    }
                    // recalculate the new fittness of this particle
                    calculateFittness(j);
                }
            }
            //after the calculation of particles in parrallel finish or in sequence
            //the herarichy starts to rearrage itself again with the new fittness
            reArrangeBreadthFirst();


        }
        /// <summary>
        /// this is helper function to get the parent of the particle returning its index
        /// </summary>
        private int getParent(int particalIndex)
        {
            int result;
            if (all_particalParentIndex[particalIndex] == -1)
            {
                result = 0;
            }
            else
            {
                result = all_particalParentIndex[particalIndex];
            }
            return result;
        }
        /// <summary>
        /// this is helper function to return the list of childern of particle returning their indexes
        /// </summary>
        private int[] getChildren(int particalIndex)
        {
            int[] result = new int[m_degree.Value];
            int counter = 0;
            for (int i = 0; i < m_ParticleSize; i++)
            {
                if (all_particalParentIndex[i] == particalIndex)
                {
                    result[counter] = i;
                    counter++;
                }
            }
            return result;
        }
        /// <summary>
        /// calculating the velocity term of the particle
        /// </summary>
        private void calculateVelocity(int particalIndex)
        {
            float rand1 = (float)R.NextDouble();
            float rand2 = (float)R.NextDouble();
            float rand3 = (float)R.NextDouble();
            
            //get the index of the parent particle
            int parentIndex = getParent(particalIndex);
            //multiple the current velocity with Inertia
            Matrix currentVelocity = (Matrix)all_velocity[particalIndex].Multiply(m_Inertia.Value);
            // calculate the diffrence between current particle and its best result found yet
            Matrix congitionPart = (Matrix)ConverttoMatrix(all_best[particalIndex]).Subtract(ConverttoMatrix(all_swarm[particalIndex]));
            congitionPart = (Matrix)congitionPart.Multiply(rand1 * m_Cof1.Value);
            // calculate the difference between the parent particle position anf the particle result
            Matrix socialPart = (Matrix)ConverttoMatrix(all_swarm[parentIndex]).Subtract(ConverttoMatrix(all_swarm[particalIndex]));
            socialPart = (Matrix)socialPart.Multiply(rand2 * m_Cof2.Value);
            
            //add to current velocity both congition part and social part
            currentVelocity = (Matrix)currentVelocity.Add(congitionPart);
            currentVelocity = (Matrix)currentVelocity.Add(socialPart);
            //clip the veolcity of the paricle and set new matrix of the particle velocity
            all_velocity[particalIndex] = clipVelocity(currentVelocity);
            all_velocity[particalIndex] = currentVelocity;
        }
        /// <summary>
        /// clip the velocity matrix
        /// </summary>
        private Matrix clipVelocity(Matrix currentVelocity)
        {
            for (int i = 0; i < currentVelocity.RowCount; i++)
            {
                for (int j = 0; j < currentVelocity.ColumnCount; j++)
                {
                    if (currentVelocity[i, j] > m_MaxVelocity)
                        currentVelocity[i, j] = m_MaxVelocity.Value;
                    if (currentVelocity[i, j] < m_MinVelocity)
                        currentVelocity[i, j] = m_MinVelocity.Value;
                }
            }
            return currentVelocity;
        }
        /// <summary>
        /// returning the level of this particle in the herarchy
        /// </summary>
        private int getParticleLevel(int particalIndex)
        {
            int counter = 0;
            for (int height = 0; height < m_height; height++)
            {
                counter += m_Particlelevel[height];
                if (counter > particalIndex)
                {
                    return height;
                }
            }
            return -1;
        }
        /// <summary>
        /// recalculate the new result from the velocity
        /// </summary>
        private void reCalculate(int particalIndex)
        {
            //get list of sites in random arragment
            int[] listofsites = GenerateRandomSolution();
            float randomGenerated;
            //get a variable of the new veolcity
            Matrix CurrentVelocity = (Matrix)all_velocity[particalIndex].Clone();
            Vector MinValues = new DenseVector(m_VariablesNumber,float.MinValue);
            ArrayList items;
            //initalize needed objects
            float sumation;
            float wheelsize;
            int currentposition;
            double result;
            // for every facility and location
            for (int i = 0; i < m_VariablesNumber; i++)
            {
                wheelsize = 0;
                sumation = 0;
                items = new ArrayList();
                //get a random first position. 
                currentposition = listofsites[i];
                randomGenerated = (float)R.NextDecimal();
                for (int j = 0; j < m_VariablesNumber; j++)
                {
                    result = 0;
                    //if this position wasnt taken allready
                    if ((CurrentVelocity[currentposition, j] != float.MinValue))
                    {
                        //calculate the sigmode value of the velocity
                        result = CurrentVelocity[currentposition, j];
                        //if the biased term is added to the system
                        if (NewbiasedCondition)
                        {
                            result += (1 / m_Matrixhurstic[currentposition, j]);
                        }
                        result = 1 / (1 + Math.Exp(-result));
                        items.Add(new retullelwheel(j, (float)result));
                    }
                }
                //sumation of the velocity values of all particles
                for (int k = 0; k < items.Count; k++)
                {
                    sumation += ((retullelwheel)items[k]).velocity;
                }
                // calculate the random probablity of this location
                for (int k = 0; k < items.Count; k++)
                {
                    ((retullelwheel)items[k]).probability = ((retullelwheel)items[k]).velocity / sumation;
                }
                if (items.Count > 1)
                {
                    for (int k = 0; k < items.Count; k++)
                    {
                        wheelsize += ((retullelwheel)items[k]).probability;
                        //if the random value pointed to this index
                        if (wheelsize >= randomGenerated)
                        {
                            // the new current position is this index
                            all_swarm[particalIndex][currentposition] = ((retullelwheel)items[k]).index;
                            // set this postion to taken
                            CurrentVelocity.SetColumn(((retullelwheel)items[k]).index, MinValues);
                            break;
                        }
                    }
                }
                else
                {
                    //the size of items are one this this is the last position
                    all_swarm[particalIndex][currentposition] = ((retullelwheel)items[0]).index;
                    CurrentVelocity.SetColumn(((retullelwheel)items[0]).index, MinValues);
                    break;
                }

            }
        }
        /// <summary>
        /// this function is to rearrage the particles in the herarchy 
        /// </summary>
        private void reArrangeBreadthFirst()
        {
            int index = -1;
            //move from the top particle
            for (int i = 0; i < m_ParticleSize - m_degree; i++)
            {
                //get the children index of the particle
                int[] childrenindex = getChildren(i);
                long max = long.MaxValue;
                //get the best fittness in the children of the particle
                for (int j = 0; j < childrenindex.Length; j++)
                {
                    if (max > all_bestfitness[childrenindex[j]])
                    {
                        index = childrenindex[j];
                        max = all_bestfitness[childrenindex[j]];
                    }
                }
                // if the best preforming children fittness better than its fittness swap particles
                if (all_bestfitness[i] < all_bestfitness[index])
                {
                    swapParticals(i, index);
                }
            }
        }
        /// <summary>
        /// this is helper function to swap particles with each other, all particle information is swapped
        /// such as the best value the best fittness the current fittness and the velocity matrix
        /// </summary>
        private void swapParticals(int i, int p)
        {
            int[] besttemp = null;
            int[] swarmtemp = null;
            Matrix velocitytemp = null;
            long fittnesstemp = 0;
            long bestfitnesstemp = 0;

            besttemp = all_best[i];
            swarmtemp = all_swarm[i];
            velocitytemp = all_velocity[i];
            fittnesstemp = all_fitness[i];
            bestfitnesstemp = all_bestfitness[i];

            all_bestfitness[i] = all_bestfitness[p];
            all_best[i] = all_best[p];
            all_swarm[i] = all_swarm[p];
            all_velocity[i] = all_velocity[p];
            all_fitness[i] = all_fitness[p];

            all_best[p] = besttemp;
            all_swarm[p] = swarmtemp;
            all_velocity[p] = velocitytemp;
            all_fitness[p] = fittnesstemp;
            all_bestfitness[p] = bestfitnesstemp;


        }
        /// <summary>
        /// this is function responsible for the local search. it search through all premutations and use the best one of them all
        /// </summary>
        private void localsearch(int particleIndex)
        {
            // starts the timer for local search
            stopwatch_LocalSearch.Start();
            
            int[] currentposition = all_swarm[particleIndex];
            int[] firstpremutation = GenerateRandomSolution();
            int[] secoundpremutation = GenerateRandomSolution();
            int index1 =0 , index2 =0, delta = 0, chosenindex1= 0, chosenindex2=0;
            int maxResult = int.MaxValue;
            
            for (int i = 0; i < m_VariablesNumber; i++)
            {
                //first index to change
                index1 = firstpremutation[i];
                for (int j = 0; j < m_VariablesNumber; j++)
                {
                    //second index to change
                    index2 = secoundpremutation[j];
                    //calculate the differnce
                    delta = calculateDiffrence(currentposition, index1, index2);
                    if (delta < maxResult && delta < 0)
                    {
                        maxResult = delta;
                        chosenindex1 = index1;
                        chosenindex2 = index2;
                    }
                }
            }
            //if better premuation found
            //get the swapped results
            if (maxResult < int.MaxValue && index1 < m_VariablesNumber && index2 < m_VariablesNumber)
            {
               all_swarm[particleIndex] = swap(chosenindex1, chosenindex2, currentposition);
            }
            stopwatch_LocalSearch.Stop();
        }
        /// <summary>
        /// this is helper function to calculate the delta of this swap
        /// </summary>
        private int calculateDiffrence(int[] currentposition, int index1, int index2)
        {
            long delta;
            delta = (m_A[index1][index1] - m_A[index2][index2]) * (m_B[currentposition[index2]][currentposition[index2]] - m_B[currentposition[index2]][currentposition[index2]])
                    + (m_A[index1][index2] - m_A[index2][index1]) * (m_B[currentposition[index2]][currentposition[index1]] - m_B[currentposition[index1]][currentposition[index2]]);
            for (int k = 0; k < m_VariablesNumber; k++)
            {
                if (k != index2 && k != index1)
                {
                    delta += (m_A[k][index1] - m_A[k][index2]) * (m_B[currentposition[k]][currentposition[index2]] - m_B[currentposition[k]][currentposition[index2]]) 
                        + (m_A[index1][k] - m_A[index2][k]) * (m_B[currentposition[index2]][currentposition[k]] - m_B[currentposition[index1]][currentposition[k]]);
                }
            }
            return (int)delta;
        }
        /// <summary>
        /// helper function in computing delta for tabu search
        /// </summary>
        public int compute_delta(int n, int[] p, int i, int j)
        {
            int d;
            int k;
            d = (m_A[i][i] - m_A[j][j]) * (m_B[p[j]][p[j]] - m_B[p[i]][p[i]]) + (m_A[i][j] - m_A[j][i]) * (m_B[p[j]][p[i]] - m_B[p[i]][p[j]]);
            for (k = 0; k < n; k = k + 1)
                if (k != i && k != j)
                    d = d + (m_A[k][i] - m_A[k][j]) * (m_B[p[k]][p[j]] - m_B[p[k]][p[i]]) + (m_A[i][k] - m_A[j][k]) * (m_B[p[j]][p[k]] - m_B[p[i]][p[k]]);
            return (d);
        }
        /// <summary>
        /// helper function for tabu search 
        /// </summary>
        public int compute_delta_part(int[] p, int[][] delta, int i, int j, int r, int s)
        {
            return (delta[i][j] + (m_A[r][i] - m_A[r][j] + m_A[s][j] - m_A[s][i]) * (m_B[p[s]][p[i]] - m_B[p[s]][p[j]] + m_B[p[r]][p[j]] - m_B[p[r]][p[i]]) + (m_A[i][r] - m_A[j][r] + m_A[j][s] - m_A[i][s]) * (m_B[p[i]][p[s]] - m_B[p[j]][p[s]] + m_B[p[j]][p[r]] - m_B[p[i]][p[r]]));
        }
        /// <summary>
        /// Tabu search function
        /// </summary>
        /// <param name="index">particle index</param>
        /// <param name="tabu_duration">tabu size</param>
        /// <param name="aspiration">Aspiration crateria</param>
        /// <param name="nr_failure">number of failures</param>
        public void tabu_search(int index, int tabu_duration, int aspiration, int nr_failure)
        {
            int[] p = new int[m_VariablesNumber];
            int[][] delta = new int[m_VariablesNumber][];
            int[][] tabu_list = new int[m_VariablesNumber][]; // tabu status
            int current_iteration; // current iteration
            int current_cost; // current sol. value
            int i; // indices
            int j;
            int k;
            int i_retained;
            int j_retained;
            int min_delta; // retained move cost
            bool autorized; // move not tabu?
            bool aspired; // move forced?
            bool already_aspired;
            long best_cost = all_fitness[index];
            int[] best_sol = new int[m_VariablesNumber];
            int failure = 0;

            for (int cont = 0; cont < m_VariablesNumber; cont++)
            {
                delta[cont] = new int[m_VariablesNumber];
                tabu_list[cont] = new int[m_VariablesNumber];
            }

            for (i = 0; i < m_VariablesNumber; i = i + 1)
            {
                p[i] = all_swarm[index][i];
                best_sol[i] = all_swarm[index][i];
            }

            current_cost = 0;
            for (i = 0; i < m_VariablesNumber; i = i + 1)
                for (j = 0; j < m_VariablesNumber; j = j + 1)
                {
                    current_cost = current_cost + m_A[i][j] * m_B[p[i]][p[j]];
                    if (i < j)
                    {
                        delta[i][j] = compute_delta(m_VariablesNumber, p, i, j);
                    }
                }
            best_cost = current_cost;

            /****************** tabu list initialization *******************/
            for (i = 0; i < m_VariablesNumber; i = i + 1)
                for (j = 0; j < m_VariablesNumber; j = j + 1)
                    tabu_list[i][j] = -(m_VariablesNumber * i + j);

            /******************** main tabu search loop ********************/
            current_iteration = 1;
            while (failure <= nr_failure)
            { //* find best move (i_retained, j_retained) *
                i_retained = int.MaxValue; // in case all moves are tabu
                j_retained = int.MaxValue;
                min_delta = int.MaxValue;
                already_aspired = false;

                for (i = 0; i < m_VariablesNumber - 1; i = i + 1)
                {
                    for (j = i + 1; j < m_VariablesNumber; j = j + 1)
                    {
                        autorized = (tabu_list[i][p[j]] < current_iteration) || (tabu_list[j][p[i]] < current_iteration);

                        aspired = (tabu_list[i][p[j]] < current_iteration - aspiration) || (tabu_list[j][p[i]] < current_iteration - aspiration) || (current_cost + delta[i][j] < best_cost);

                        if ((aspired && !already_aspired) || (aspired && already_aspired && (delta[i][j] < min_delta)) ||
                           (!aspired && !already_aspired && (delta[i][j] < min_delta) && autorized))
                        {
                            i_retained = i; j_retained = j;
                            min_delta = delta[i][j];
                            if (aspired) { already_aspired = true; };
                        }
                    }
                }

                if (i_retained != int.MaxValue)
                { //* transpose elements in pos. i_retained and j_retained *

                    p = swap(i_retained, j_retained, p);
                    /* update solution value*/
                    current_cost = current_cost + delta[i_retained][j_retained];
                    /* forbid reverse move for a random number of iterations*/
                    tabu_list[i_retained][p[j_retained]] = current_iteration + (int)(Math.Pow(R.NextDouble(), 3) * tabu_duration);
                    tabu_list[j_retained][p[i_retained]] = current_iteration + (int)(Math.Pow(R.NextDouble(), 3) * tabu_duration);

                    /* best solution improved ?*/
                    if (current_cost < best_cost)
                    {
                        best_cost = current_cost;
                        for (k = 0; k < m_VariablesNumber; k = k + 1)
                            best_sol[k] = p[k];
                    }
                    else
                    {
                        failure++;
                    }

                    /* update matrix of the move costs*/
                    for (i = 0; i < m_VariablesNumber - 1; i = i + 1)
                    {
                        for (j = i + 1; j < m_VariablesNumber; j = j + 1)
                        {
                            if (i != i_retained && i != j_retained && j != i_retained && j != j_retained)
                            {
                                delta[i][j] = compute_delta_part(p, delta, i, j, i_retained, j_retained);
                            }
                            else
                            {
                                delta[i][j] = compute_delta(m_VariablesNumber, p, i, j);
                            }
                        }
                    }
                }

                current_iteration = current_iteration + 1;

            }

            for (i = 0; i < m_VariablesNumber; i = i + 1)
            {
                all_swarm[index][i] = best_sol[i];
            }
            all_fitness[index] = best_cost;
        }
        /// <summary>
        /// progress report as final report of the system
        /// </summary>
        private void reportProgree()
        {
            DirectoryInfo directory = new DirectoryInfo(DirectoryName);
            if (!directory.Exists)
            {
                directory.Create();
            }
            string currentname = "HPSO";
            if (fulllocal)
                currentname += "FullyLocalSearch";
            if (singlelocal)
                currentname += "SingleLocalSearch";
            if (NewbiasedCondition)
                currentname += "NewbiasedTerm";
            if (simulatedAnnealing)
                currentname += "badswitch";
            if (Fulltabou)
                currentname += "FullTabou";
            if (singletabou)
                currentname += "singleTabou";
            DirectoryInfo[] listofdirectories = directory.GetDirectories();
            bool Exist = false;
            for (int i = 0; i < listofdirectories.Length; i++)
            {
                if (listofdirectories[i].Name == currentname)
                {
                    Exist = true;
                }
            }
            if (!Exist)
            {
                directory = directory.CreateSubdirectory(currentname);
            }
            else
            {
                directory = new DirectoryInfo(DirectoryName + "\\" + currentname);
            }
            string[] arr = m_FileNameRead.Split('\\');
            string filename = "";
            filename += Guid.NewGuid().ToString();
            filename += arr[arr.Length - 1];
            StreamWriter sw = new StreamWriter(directory.FullName +  "\\" + filename);
            sw.WriteLine(string.Format("Swarm Size {0}", m_ParticleSize));
            int functionCounterIterationCounter = 0;
            sw.WriteLine();
            sw.WriteLine("Summary");
            for (int i = 0; i < iterationNumber.Count; i++)
            {
                functionCounterIterationCounter = System_functionCounter[i] * iterationNumber[i];
                sw.WriteLine(iterationNumber[i].ToString() + "\t" + bestSystemfittness[i].ToString() + "\t" + System_functionCounter[i].ToString() + "\t" + functionCounterIterationCounter.ToString());
                sw.WriteLine("Result");
                for (int j = 0; j < m_VariablesNumber; j++)
                {
                    sw.Write(System_myBestResult[i][j]);
                    sw.Write(",");
                }
                sw.WriteLine();
            }
            sw.WriteLine();
            sw.WriteLine(string.Format(String.Format("Fittness : " + m_bestSystemFitness.Value)));
            sw.WriteLine(string.Format(String.Format("intializeSystem : " + intializeSystem)));
            sw.WriteLine(string.Format(String.Format("LocalSearch : " + localSearch)));
            sw.WriteLine(string.Format(String.Format("HuristicTerm : " + huristicTerm)));
            sw.WriteLine(string.Format(String.Format("Runs : " + runs)));
            sw.WriteLine(string.Format(String.Format("Tabu Search : " + tabuSearch)));
            sw.WriteLine(string.Format(String.Format("TotalRun : " + totalRun)));
            sw.WriteLine();
            sw.Close();
        }
        /// <summary>
        /// dispose and clean the objects behind it.
        /// </summary>
        public void Dispose()
        {
            functionCounter = null;
            System_functionCounter = null;
            all_best = null;
            all_swarm = null;
            all_particalParentIndex = null;
            all_velocity = null;
            all_fitness = null;
            all_bestfitness = null;
            m_Cof1 = null;
            m_Cof2 = null;
            m_MaxVelocity = null;
            m_MinVelocity = null;
            m_height = null;
            m_degree = null;
            m_bestSystemFitness = null;
            bestSystemfittness = null;
            iterationNumber = null;
            m_ParticleSize = null;
            m_Particlelevel = null;
            GC.Collect();
            GC.WaitForPendingFinalizers();
            GC.Collect();
        }
    }
    /// <summary>
    /// helper class for the recalculation of the result
    /// </summary>
    public class retullelwheel:IComparable
    {
        public int index;
        public float velocity;
        public float probability;

        public retullelwheel(int m_index, float m_velocity)
        {
            index = m_index;
            velocity = m_velocity;
        }
        public int CompareTo(object obj)
        {
            retullelwheel objstruct = (retullelwheel)obj;
            if (objstruct.probability > probability)
                return 1;
            if (objstruct.probability < probability)
                return -1;
            return 0;
        }
    }
}
